This is a documentation for analyses of scATAC-seq data, generated from rat metrial gland tissues on gestational day (GD) 15.5 and 19.5.

GD15.5:

library(Seurat)
library(Signac)
library(ggplot2)
library(GenomicRanges)
set.seed(1234)


load("/work/LAS/geetu-lab/hhvu/project3_scATAC/rnor6.rda")
colnames(rnor6)[1] <- "gene_id"
rnor6.ranges <- makeGRangesFromDataFrame(rnor6, seqnames.field = "chromosome_name",
                                         start.field = "start_position",
                                         end.field = "end_position",
                                         strand.field = "strand",
                                         keep.extra.columns = TRUE)

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/scRNA-seqFiles/gd15.5_res0.8.rda")
data <- RenameIdents(data, `0` = "Other cells", `1` = "Other cells",
                     `2` = "Other cells", `3` = "Other cells", `4` = "Other cells",
                     `5` = "Other cells", `6` = "Natural killer cells",
                     `7` = "Macrophage cells", `8` = "Macrophage cells",
                     `9` = "Natural killer cells", `10` = "Endothelial cells",
                     `11` = "Macrophage cells", `12` = "Other cells",
                     `13` = "Endothelial cells", `14` = "Smooth muscle cells",
                     `15` = "Endothelial cells", `16` = "Other cells",
                     `17` = "Smooth muscle cells", `18` = "Other cells",
                     `19` = "Natural killer cells", `20` = "Other cells",
                     `21` = "Other cells", `22` = "Other cells",
                     `23` = "Invasive trophoblast cells", `24` = "Other cells",
                     `25` = "Other cells", `26` = "Smooth muscle cells",
                     `27` = "Other cells", `28` = "Other cells")
Idents(data) <- factor(Idents(data), levels = c("Other cells", "Macrophage cells", "Smooth muscle cells",
                                                "Endothelial cells", "Natural killer cells", "Invasive trophoblast cells"))
rats.RNA <- data

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/2_mergedSamples/2_MACS2Peaks/GD15.5/15.5_1-2-3-merged.rda")
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/rnor6.ranges.rda")

rats <- unintegrated
DefaultAssay(rats) <- "combined"
mcols(rnor6.ranges, level="within")[, "gene_id"] <- mcols(rnor6.ranges, level="within")[, "ensembl_gene_id"]
Annotation(rats) <- rnor6.ranges
gene.activities <- GeneActivity(rats)
rats[['RNA']] <- CreateAssayObject(counts = gene.activities[, colnames(rats)])
#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')
rats <- NormalizeData(object = rats, assay = 'RNA', normalization.method = 'LogNormalize', scale.factor = median(rats$nCount_RNA))
DefaultAssay(rats) <- 'RNA'

transfer.anchors <- FindTransferAnchors(reference = rats.RNA, query = rats, reduction = 'cca')
predicted.labels <- TransferData(anchorset = transfer.anchors, refdata = Idents(rats.RNA), weight.reduction = rats[['lsi']], dims=2:20)

rats <- AddMetaData(object = rats, metadata= predicted.labels)

plot1 <- DimPlot(object = rats.RNA, cols = c("grey70", "orange3", "#718748", "#183D61", "#79704C", "red3")) + ggtitle('GD15.5 scRNA-seq')
plot2 <- DimPlot(object=rats, group.by = 'predicted.id', label = T, repel = T) + ggtitle('GD15.5 scATAC-seq')

plot1 + plot2


#save(rats, file="/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/gd15.5-ATACwithRNAlables.rda")

table(rats$predicted.id)
#> 
#>          Endothelial cells Invasive trophoblast cells 
#>                       3250                        128 
#>           Macrophage cells       Natural killer cells 
#>                       2682                       1443 
#>                Other cells        Smooth muscle cells 
#>                      16700                       1118

Correlation check:

library(Seurat)
library(Signac)
library(Seurat)
library(S4Vectors)
library(GenomicRanges)
library(patchwork)
library(biomaRt)
library(ggplot2)
library(ggpubr)
#> Warning: package 'ggpubr' was built under R version 4.0.3
library(BSgenome.Rnorvegicus.Ensembl.rn6)
set.seed(1234)

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd15.5-ATACwithRNAlables.rda")
Idents(rats) <- rats@meta.data$predicted.id # change from old identities to those predicted by scRNA-seq
DefaultAssay(rats) <- 'MACSpeaks'
avgATAC <- AverageExpression(rats, assays = "RNA")
avgATAC <- as.data.frame(avgATAC)

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/scRNA-seqFiles/gd15.5_res0.8.rda")
avgRNA <- AverageExpression(data, assays = "RNA")
avgRNA <- as.data.frame(avgRNA)

Trophoblast:

avgRNA1 <- avgRNA$RNA.23
names(avgRNA1) <- rownames(avgRNA)

avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Invasive.trophoblast.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]

cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#> 
#>  Spearman's rank correlation rho
#> 
#> data:  avgRNA1 and avgATAC1
#> S = 4.4498e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.4378357

df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-15.5-tbcCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Invasive Trophoblast Cells at gd 15.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")

#dev.off()

Endothelial:

avgRNA1 <- rowMeans(avgRNA[, c("RNA.10", "RNA.13", "RNA.15")])
names(avgRNA1) <- rownames(avgRNA)

avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Endothelial.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]

cor.test(log10(avgRNA1+10^(-5)), log10(avgATAC1+10^(-5)), method = "spearman", exact = FALSE)
#> 
#>  Spearman's rank correlation rho
#> 
#> data:  log10(avgRNA1 + 10^(-5)) and log10(avgATAC1 + 10^(-5))
#> S = 3.9856e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.4964745

df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-15.5-endoCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=FALSE) + ggtitle("Correlation in Endothelial Cells at gd 15.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")

#dev.off()

Macrophage:

avgRNA1 <- rowMeans(avgRNA[, c("RNA.7", "RNA.8", "RNA.11")])
names(avgRNA1) <- rownames(avgRNA)

avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Macrophage.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]

cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#> 
#>  Spearman's rank correlation rho
#> 
#> data:  avgRNA1 and avgATAC1
#> S = 3.6982e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.5327817

df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-15.5-macroCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Macrophage Cells at gd 15.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")

#dev.off()

NK:

avgRNA1 <- rowMeans(avgRNA[, c("RNA.6", "RNA.9", "RNA.19")])
names(avgRNA1) <- rownames(avgRNA)

avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Natural.killer.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]

cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#> 
#>  Spearman's rank correlation rho
#> 
#> data:  avgRNA1 and avgATAC1
#> S = 4.0055e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.4939641

df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-15.5-nkCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Natural Killer Cells at gd 15.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")

#dev.off()

Smooth muscle:

avgRNA1 <- rowMeans(avgRNA[, c("RNA.14", "RNA.17", "RNA.26")])
names(avgRNA1) <- rownames(avgRNA)

avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Smooth.muscle.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]

cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#> 
#>  Spearman's rank correlation rho
#> 
#> data:  avgRNA1 and avgATAC1
#> S = 4.382e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.4464026

df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-15.5-smCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Smooth Muscle Cells at gd 15.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")

#dev.off()

GD19.5:

set.seed(1234)

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/scRNA-seqFiles/gd19.5-4-5-6-7_res0.8.rda")
data <- data2
data <- RenameIdents(data, `0` = "Other cells", `1` = "Macrophage cells",
                     `2` = "Other cells", `3` = "Other cells", `4` = "Other cells",
                     `5` = "Other cells", `6` = "Invasive trophoblast cells",
                     `7` = "Macrophage cells", `8` = "Smooth muscle cells",
                     `9` = "Natural killer cells", `10` = "Endothelial cells",
                     `11` = "Macrophage cells", `12` = "Other cells",
                     `13` = "Other cells", `14` = "Macrophage cells",
                     `15` = "Other cells", `16` = "Other cells",
                     `17` = "Macrophage cells", `18` = "Macrophage cells",
                     `19` = "Other cells", `20` = "Endothelial cells",
                     `21` = "Other cells", `22` = "Other cells",
                     `23` = "Other cells", `24` = "Other cells",
                     `25` = "Other cells", `26` = "Smooth muscle cells")
Idents(data) <- factor(Idents(data), levels = c("Other cells", "Macrophage cells", "Smooth muscle cells",
                                                "Endothelial cells", "Natural killer cells", "Invasive trophoblast cells"))
rats.RNA <- data

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/2_mergedSamples/2_MACS2Peaks/GD19.5/19.5_4-5-6-merged.rda")
rats <- unintegrated
DefaultAssay(rats) <- "combined"
mcols(rnor6.ranges, level="within")[, "gene_id"] <- mcols(rnor6.ranges, level="within")[, "ensembl_gene_id"]
Annotation(rats) <- rnor6.ranges
gene.activities <- GeneActivity(rats)
rats[['RNA']] <- CreateAssayObject(counts = gene.activities[, colnames(rats)])
#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')
rats <- NormalizeData(object = rats, assay = 'RNA', normalization.method = 'LogNormalize', scale.factor = median(rats$nCount_RNA))
DefaultAssay(rats) <- 'RNA'

transfer.anchors <- FindTransferAnchors(reference = rats.RNA, query = rats, reduction = 'cca')
#> Warning in RunCCA.Seurat(object1 = reference, object2 = query, features =
#> features, : Running CCA on different assays
predicted.labels <- TransferData(anchorset = transfer.anchors, refdata = Idents(rats.RNA), weight.reduction = rats[['lsi']], dims=2:10)

rats <- AddMetaData(object = rats, metadata = predicted.labels)

plot1 <- DimPlot(object = rats.RNA, cols = c("grey70", "orange3", "#718748", "#183D61", "#79704C", "red3")) + ggtitle('GD19.5 scRNA-seq')
plot2 <- DimPlot(object=rats, group.by = 'predicted.id', label = T, repel = T) + ggtitle('GD19.5 scATAC-seq')

plot1 + plot2


#save(rats, file="/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/gd19.5-ATACwithRNAlables.rda")

table(rats$predicted.id)
#> 
#>          Endothelial cells Invasive trophoblast cells 
#>                        894                        684 
#>           Macrophage cells       Natural killer cells 
#>                       1265                        223 
#>                Other cells        Smooth muscle cells 
#>                      11125                        197

Correlation check:

set.seed(1234)

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-ATACwithRNAlables.rda")
Idents(rats) <- rats@meta.data$predicted.id # change from old identities to those predicted by scRNA-seq
DefaultAssay(rats) <- 'MACSpeaks'
avgATAC <- AverageExpression(rats, assays = "RNA")
avgATAC <- as.data.frame(avgATAC)

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/scRNA-seqFiles/gd19.5-4-5-6-7_res0.8.rda")
avgRNA <- AverageExpression(data2, assays = "RNA")
avgRNA <- as.data.frame(avgRNA)

Trophoblast:

avgRNA1 <- avgRNA$RNA.6
names(avgRNA1) <- rownames(avgRNA)

avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Invasive.trophoblast.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]

cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#> 
#>  Spearman's rank correlation rho
#> 
#> data:  avgRNA1 and avgATAC1
#> S = 4.109e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.5258185

df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-19.5-tbcCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho") + ggtitle("Correlation in Invasive Trophoblast Cells at gd 19.5")

#dev.off()

Endothelial:

avgRNA1 <- rowMeans(avgRNA[, c("RNA.10", "RNA.20")])
names(avgRNA1) <- rownames(avgRNA)

avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Endothelial.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]

cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#> 
#>  Spearman's rank correlation rho
#> 
#> data:  avgRNA1 and avgATAC1
#> S = 3.9959e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.5388751

df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-19.5-endoCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Endothelial Cells at gd 19.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")

#dev.off()

Macrophage:

avgRNA1 <- rowMeans(avgRNA[, c("RNA.1", "RNA.7", "RNA.11", "RNA.17", "RNA.18")])
names(avgRNA1) <- rownames(avgRNA)

avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Macrophage.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]

cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#> 
#>  Spearman's rank correlation rho
#> 
#> data:  avgRNA1 and avgATAC1
#> S = 4.0517e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>      rho 
#> 0.532435

df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-19.5-macroCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Macrophage Cells at gd 19.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")

#dev.off()

NK:

avgRNA1 <- avgRNA[, c("RNA.9")]
names(avgRNA1) <- rownames(avgRNA)

avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Natural.killer.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]

cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#> 
#>  Spearman's rank correlation rho
#> 
#> data:  avgRNA1 and avgATAC1
#> S = 4.0328e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.5346098

df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-19.5-nkCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Natural Killer Cells at gd 19.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")

#dev.off()

Smooth muscle:

avgRNA1 <- rowMeans(avgRNA[, c("RNA.8", "RNA.26")])
names(avgRNA1) <- rownames(avgRNA)

avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Smooth.muscle.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]

cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#> 
#>  Spearman's rank correlation rho
#> 
#> data:  avgRNA1 and avgATAC1
#> S = 4.438e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.4878489

df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-19.5-smCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Smooth Muscle Cells at gd 19.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")

#dev.off()
sessionInfo()
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#>  [1] BSgenome.Rnorvegicus.Ensembl.rn6_0.0.1
#>  [2] BSgenome_1.56.0                       
#>  [3] rtracklayer_1.48.0                    
#>  [4] Biostrings_2.56.0                     
#>  [5] XVector_0.28.0                        
#>  [6] ggpubr_0.4.0                          
#>  [7] biomaRt_2.44.4                        
#>  [8] patchwork_1.0.1                       
#>  [9] GenomicRanges_1.40.0                  
#> [10] GenomeInfoDb_1.24.2                   
#> [11] IRanges_2.22.2                        
#> [12] S4Vectors_0.26.1                      
#> [13] BiocGenerics_0.34.0                   
#> [14] ggplot2_3.3.6                         
#> [15] Signac_1.5.0                          
#> [16] SeuratObject_4.0.4                    
#> [17] Seurat_4.1.0                          
#> 
#> loaded via a namespace (and not attached):
#>  [1] utf8_1.2.2             reticulate_1.16        tidyselect_1.1.2      
#>  [4] RSQLite_2.2.1          AnnotationDbi_1.50.3   htmlwidgets_1.5.2     
#>  [7] grid_4.0.2             docopt_0.7.1           BiocParallel_1.22.0   
#> [10] Rtsne_0.15             munsell_0.5.0          codetools_0.2-16      
#> [13] ica_1.0-2              future_1.19.1          miniUI_0.1.1.1        
#> [16] withr_2.5.0            colorspace_2.0-3       Biobase_2.48.0        
#> [19] highr_0.9              knitr_1.39             rstudioapi_0.13       
#> [22] ROCR_1.0-11            ggsignif_0.6.0         tensor_1.5            
#> [25] listenv_0.8.0          labeling_0.4.2         slam_0.1-47           
#> [28] GenomeInfoDbData_1.2.3 polyclip_1.10-0        bit64_4.0.5           
#> [31] farver_2.1.0           vctrs_0.4.1            generics_0.1.2        
#> [34] xfun_0.31              BiocFileCache_1.12.1   lsa_0.73.2            
#> [37] ggseqlogo_0.1          R6_2.5.1               DelayedArray_0.14.1   
#> [40] bitops_1.0-6           spatstat.utils_2.2-0   assertthat_0.2.1      
#> [43] promises_1.1.1         scales_1.2.0           gtable_0.3.0          
#> [46] globals_0.13.0         goftest_1.2-2          rlang_1.0.3           
#> [49] RcppRoll_0.3.0         splines_4.0.2          rstatix_0.6.0         
#> [52] lazyeval_0.2.2         spatstat.geom_2.3-0    broom_0.8.0           
#> [55] yaml_2.3.5             reshape2_1.4.4         abind_1.4-5           
#> [58] backports_1.4.1        httpuv_1.5.4           tools_4.0.2           
#> [61] ellipsis_0.3.2         spatstat.core_2.3-1    jquerylib_0.1.4       
#> [64] RColorBrewer_1.1-3     ggridges_0.5.2         Rcpp_1.0.8            
#> [67] plyr_1.8.6             progress_1.2.2         zlibbioc_1.34.0       
#> [70] purrr_0.3.4            RCurl_1.98-1.2         prettyunits_1.1.1     
#> [73] rpart_4.1-15           openssl_2.0.2          deldir_1.0-6          
#>  [ reached getOption("max.print") -- omitted 88 entries ]